% This script calculates PNJunction of Graphene 
% Physical Constant Setup (Hartree Atomic Unit)
eps = 1;               % Relative permittivity
eps_0 = 1/4/pi;   % Vacuum permittivity 
e = 1; % Electronic Charge
hbar = 1;
Vh = 27.2113850560;
kB = 8.6173324E-5/27.21138386;
vf = 1E6/(2.187691263373E6);

% Energy Grid
Ne = 1500; nec = (Ne)/2;
Emax = 0.05;
Emin = -0.05;
he = (Emax - Emin)/(Ne-1);
E = linspace(Emin, Emax, Ne);


% DOS1 (In Hartree Atomic Unit)
D1 = (2/pi/hbar^2/vf^2)*abs(E);
Ev1 = 0; Ec1 = 0;

% DOS2 (In Hartree Atomic Unit)
D2 = (2/pi/hbar^2/vf^2)*abs(E);
Ev2 = 0; Ec2 = 0;

% Spatial Grid (In Hartree Atomic Unit)
L = 4000;
Nx = 4096;
ncx = (Nx)/2;
x = linspace(-L/2,L/2,Nx);

% Charge Setup
f1 = -1E-4;
f2 = +5E-4;
Q1 = f1 + sum(D1(1:round((Ev1-Emin)/he)))*he;
Q2 = f2 + sum(D2(1:round((Ev2-Emin)/he)))*he;

%% Initial Calculation
[Va, Ra, Err] = PNJunction(Q1, Q2, [E;D1], [E;D2], 100*kB, 5E-3, Nx, L);
filename = sprintf('pnGrNew(f1=%1.1E,f2=%1.1E,L=%i,N=%i)',f1,f2,L,Nx);
dlmwrite([filename,'V.txt'],[x',Va']);
dlmwrite([filename,'R.txt'],[x',Ra']);
dlmwrite([filename,'Err.txt'],Err);

%% Continued Calculation
filename = sprintf('pnGrNew(f1=%1.1E,f2=%1.1E,L=%i,N=%i)',f1,f2,L,Nx);
Vin = dlmread([filename,'V.txt']);
[Va, Ra, Err] = PNJunction(Q1, Q2, [E;D1], [E;D2], 10*kB, 5E-3, Nx, L, Vin(:,2)');
dlmwrite([filename,'V.txt'],[x',Va']);
dlmwrite([filename,'R.txt'],[x',Ra']);
dlmwrite([filename,'Err.txt'],Err);

%% pngraphene Calculation
[V0,R0,Err0] = PNGraphene(f2,1E-3,Nx,L);
filename = sprintf('pnGr(f=%1.1E,L=%i,N=%i)',f2,L,Nx);
dlmwrite([filename,'V.txt'],[x',V0']);
dlmwrite([filename,'R.txt'],[x',R0']);
dlmwrite([filename,'Err.txt'],Err0);


%% Comparison
filename0 = sprintf('pnGr(f=%1.1E,L=%i,N=%i)',f2,L,Nx);
filename1 = sprintf('pnGrNew(f=%1.1E,L=%i,N=%i)',f2,L,Nx);
V0 = dlmread([filename0,'V.txt']);
R0 = dlmread([filename0,'R.txt']);
V1 = dlmread([filename1,'V.txt']);
R1 = dlmread([filename1,'R.txt']);

figure; 
subplot(2,1,1); hold on;
plot(V0(:,1),V0(:,2),'-k');
plot(V1(:,1),V1(:,2),'-r');

subplot(2,1,2); hold on;
plot(R0(:,1),1./R0(:,2),'-k');
plot(R1(:,1),1./R1(:,2),'-r');
xlim([1,500]);

%% Plot Asymmetric Result
filename = sprintf('pnGrNew(f1=%1.1E,f2=%1.1E,L=%i,N=%i)',f1,f2,L,Nx);
V = dlmread([filename,'V.txt']);
R = dlmread([filename,'R.txt']);

figure;
subplot(2,1,1); hold on;
plot(V(:,1),V(:,2),'-k');
xlim([-100,100]);

subplot(2,1,2); hold on;
x = R(:,1);
R_metal = (0.0337)/(2*pi^2)./x;
plot(x,R_metal,'-r');
plot(R(:,1),R(:,2),'-k');
ylim([-4E-4,8E-4]);
xlim([-50,50]);


